#! /usr/bin/python
# -*- coding: utf-8 -*-
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#
#  Create the data files to be used to test the program qview.
#_________________________________________________________
# USAGE:
#    ./makeData.py
#_________________________________________________________
# A.Z. - 10/07 => Creation
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
import sys
import numpy


# Time interval (s)
t = numpy.arange(  0.0, 1000.0, 1.0, dtype=float )
Nt= len(t)


# Light Curves (I'm not so creative as you can see!)
p = [] # Periods
p.append(100.0) # s
p.append(300.0) # s
p.append(500.0) # s
p.append(700.0) # s
a = []  # Amplitudes
a.append( 500.0) # counts/s
a.append(1000.0) # counts/s
a.append(1500.0) # counts/s
a.append(2000.0) # counts/s
sn= []  # Signal to noise
sn.append( 3.0)
sn.append( 4.0)
sn.append( 5.0)
sn.append( 6.0)

# Light Curves
N = len(p)
y = []
for i in range(0,N):
  w = 2.0*numpy.pi/p[i]  # Angular Frequency
  y.append(  10000.0 + a[i]*numpy.sin(w*t) )
  y[i] = y[i] + (a[i]/sn[i])*numpy.random.normal(0.0, 1.0, Nt)




# Save light curves
header ="# COLUMN  1: Time (Julian Date like)\n"
header+="# COLUMN >2: Counts (s^-1)\n"
out = file("lc.dat",'w')
out.write(header)
for i in range(0,Nt):
  # Fake Julian date. Just to have the formated light curve as in a real
  # observation
  tj = 2454552.0 + t[i]/(24.0*60.0*60.0)

  out.write( "%.10f  " % tj )
  for j in range(0,N):
    out.write( "%.6f  " % y[j][i] )
  out.write("\n")
out.close()


from pylab import *
for i in range(0,N):
  plot(t,y[i],'o')
show()
